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The  Quantum  Normal  Form  Approach  to  Reactive  Scattering:  The  Cumulative 
Reaction  Probability  for  Collinear  Exchange  Reactions 


Arseni  Goussev^,  Roman  Schubert^,  Holger  Waalkens^’^,  and  Stephen  Wiggins^ 
^School  of  Mathematics,  University  of  Bristol,  University  Walk,  Bristol  BS8  ITW,  UK 
^Department  of  Mathematics  and  Computing  Science, 

University  of  Groningen,  PO  Box  407,  9700  AK  Groningen,  The  Netherlands 

(Dated:  August  5,  2009) 

The  quantum  normal  form  approach  to  quantum  transition  state  theory  is  used  to  compute  the 
cumulative  reaction  probability  for  collinear  exchange  reactions.  It  is  shown  that  for  heavy  atom 
systems  like  the  nitrogen  exchange  reaction  the  quantum  normal  form  approach  gives  excellent 
results  and  has  major  computational  benefits  over  full  reactive  scattering  approaches.  For  light 
atom  systems  like  the  hydrogen  exchange  reaction  however  the  quantum  normal  approach  is  shown 
to  give  only  poor  results.  This  failure  is  attributed  to  the  importance  of  tunnelling  trajectories  in 
light  atom  reactions  that  are  not  captured  by  the  quantum  normal  form  as  indicated  by  the  only 
very  slow  convergence  of  the  quantum  normal  form  for  such  systems. 

PACS  numbers:  34.10.+X,  34.50.Lf,  05.,  02.70.-c,  02.,03.65.Xp,82.20.-w,82.20.Db,82.20.Ej,82.30.Hk 


I.  INTRODUCTION 

The  classical  mechanical  picture  of  a  chemical  reaction 
as  a  scattering  problem  across  a  saddle  point  of  the  Born- 
Oppenheimer  potential  energy  surface  in  configuration 
space  has  proven  to  be  a  fruitful  way  of  visualizing  and 
thinking  about  chemical  reactions  since  the  1930’s,  when 
Eyring,  Polanyi,  and  Wigner  developed  transition  state 
theory  (TST).  TST  provides  the  framework  for  comput¬ 
ing,  using  classical  mechanics,  many  of  the  physically  im¬ 
portant  quantities  for  describing  such  chemical  reactions. 
The  fundamental  geometrical  object  in  TST  is  a  divid¬ 
ing  surface  that  divides  the  energy  surface  into  a  reactant 
and  a  product  component.  With  such  a  dividing  surface 
in  hand,  one  can  then  compute  the  reaction  rate  from  the 
directional  phase  space  flux  through  this  surface.  In  or¬ 
der  not  to  overestimate  the  rate  the  dividing  surface  must 
not  be  recrossed  by  reactive  trajectories,  i.e.  the  dividing 
surface  should  have  the  “no  re-crossing”  property.  In  the 
70’s  Pechukas,  Poliak  and  others  [1,  2]  showed  that  for 
two  degrees  of  freedom  such  a  dividing  surface  can  be  con¬ 
structed  from  a  periodic  orbit  (the  so  called  periodic  or¬ 
bit  dividing  surface).  Recently  it  has  been  shown  that  for 
more  than  two  degrees-of-freedom  a  dividing  surface  that 
is  free  of  recrossings  can  be  built  from  a  normally  hyper¬ 
bolic  invariant  manifold  (NHIM)  [3].  The  dividing  sur¬ 
face  and  the  NHIM  can  be  directly  constructed  from  an 
algorithm  based  on  a  Poincare-Birkhoff  normal  form  pro¬ 
cedure  [4]  which  also  gives  an  expression  for  the  flux  [5]. 
The  classical  phase  space  transition  state  theory,  based 
on  Poincare-Birkhoff  normal  form  theory,  naturally  leads 
to  a  quantum  version  of  transition  state  theory,  based  on 
a  quantum  normal  form.  Since  the  normal  form  is  valid  in 
a  neighborhood  in  energy  both  above  and  below  the  sad¬ 
dle  point,  it  includes  the  quantum  effect  of  tunneling  in 
the  region  near  the  saddle.  Moreover,  it  does  not  require 
a  full  quantum  simulation  in  a  neighborhood  of  the  TST 
dividing  surface  ([6,  7])  in  order  to  compute  important 
quantities  associated  with  the  reaction.  This  is  signih- 


cant  since  much  effort  has  been  devoted  to  developing  a 
quantum  version  of  transition  state  theory  whose  imple¬ 
mentation  remains  feasible  for  multi-dimensional  systems 
(see  the  flux-flux  autocorrelation  function  formalism  by 
Miller  and  coworkers  [8]).  However,  in  [9]  Miller  stated 
that  “-the  conclusion  of  it  all  is  that  there  is  no  uniquely 
well  defined  quantum  version  of  TST  in  the  sense  that 
there  is  in  classical  mechanics.  This  is  because  tunnel¬ 
ing  along  the  reaction  coordinate  necessarily  requires  one 
to  solve  the  (quantum)  dynamics  for  some  hnite  region 
about  the  TS  dividing  surface,  and  if  one  does  this  quan¬ 
tum  mechanically  there  is  no  ‘theory’  left,  i.e.,  one  has 
a  full  dimensional  quantum  dynamics  treatment  that  is 
ipso  facto  exact,  a  quantum  simulation.”  Nevertheless, 
our  approach  based  on  the  quantum  normal  form  leads 
to  a  quantum  version  of  transition  state  theory  that  in¬ 
cludes  tunneling  near  the  saddle  and  does  not  require  a 
full  quantum  simulation  in  a  neighborhood  of  the  TST  di¬ 
viding  surface.  Moreover,  our  computation  of  the  cumu¬ 
lative  reaction  probability  can  be  viewed  as  the  quantum 
mechanical  flux  through  a  (classically  recrossing  free)  di¬ 
viding  surface,  which  includes  tunneling.  The  quantum 
normal  form  gives  a  local  decoupling  of  the  quantum  dy¬ 
namics  to  any  desired  order  in  ft  which  is  the  key  issue 
here,  i.e.  locally,  we  have  a  decoupling  of  the  scattering 
states  into  forward/backward  reactive  and  non-reactive, 
and  for  these  states  we  know  the  transmission  probabili¬ 
ties  analytically.  Therefore  we  do  not  have  to  ‘simulate’ 
the  quantum  dynamics.  Hence,  it  sidesteps  the  issues 
and  concerns  expressed  by  Miller. 

In  this  paper  we  illustrate  the  utility  of  the  quantum 
normal  form  approach  to  quantum  transition  state  theory 
by  considering  the  computation  and  behavior  of  the  bi- 
molecular  cumulative  reaction  probability  (CRP)  Af{E), 
dehned  as  [10,  11] 


A7(A)  =  ^  \Sn.,n,{E)f,  (1) 

nr,np 
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where  S{E)  is  the  reactive  scattering  matrix  evaluated  at 
energy  E,  and  nj.  (np)  are  the  quantum  numbers  describ¬ 
ing  the  asymptotic  channel  of  incoming  reactants  (outgo¬ 
ing  products).  The  CRP  is  a  fundamental  quantity  that 
characterizes  the  reaction  rate:  the  microcanonical  and 
canonical  rate  constants  can  be  determined  from  Af{E) 
by  means  of  simple  relations  [10]. 

This  paper  is  outlined  as  follows.  In  Section  II  we 
outline  the  theoretical  and  computational  aspects  of  the 
quantum  normal  form  theory.  We  emphasize  the  struc¬ 
tural  features  that  allow  for  the  treatment  of  high  di¬ 
mensional  quantum  problems  and  show  how  the  quan¬ 
tum  normal  form  leads  to  a  “simple”  expression  for  the 
CRP.  In  Section  III  we  apply  the  quantum  normal  form 
approach  to  the  computation  of  the  CRP  for  the  collinear 
hydrogen  and  nitrogen  exchange  reactions.  These  quan¬ 
tities  are  compared  to  the  “exact”  answer  obtained  from 
a  reactive  quantum  scattering  calculation.  In  Section  IV 
we  discuss  some  aspects  related  to  the  convergence  prop¬ 
erties  of  the  quantum  normal  form,  and  in  Section  V  we 
summarize  our  results  and  offer  some  directions  for  fur¬ 
ther  investigations. 


II.  QUANTUM  NORMAL  FORM  THEORY 

In  this  section  we  present  central  aspects  of  the  quan¬ 
tum  normal  form  (QNF)  theory;  for  rigorous  mathemat¬ 
ical  statements,  proofs  and  further  details  we  refer  to 
Ref.  [6,  7]. 

We  begin  by  considering  a  quantum  Hamilton  operator 
E[  which  we  assume  to  be  obtained  from  the  Weyl  quan¬ 
tization  of  a  classical  Hamilton  function  iJ(p,  q).  Here 
q  =  and  p  =  (j)i,p2,  ■ .  ■  ,Pd)  denote  the 

canonical  coordinates  and  momenta,  respectively,  of  a 
Hamiltonian  system  with  d  degrees  of  freedom.  Through¬ 
out  this  paper  we  will  use  atomic  units,  so  that  q  and  p 
are  dimensionless.  We  will  denote  the  corresponding  op¬ 
erators  by  q  =  iqi,h,  ■■■,qd)  and  p  =  ipi,P2,  ■  ■  ■  ,Pd)- 
In  the  coordinate  representation  their  components  corre¬ 
spond  to  multiplication  by  Qj  and  the  differential  opera¬ 
tors  pj  =  —ihesd/dqj.  Here  hes  is  a  dimensionless  pa¬ 
rameter  which  corresponds  to  a  scaled,  effective  Planck’s 
constant.  For  molecular  reactions  described  in  the  Born- 
Oppenheimer  approximation,  occurs  naturally  as  the 
ratio  of  the  electronic  mass  and  the  reduced  mass  of  the 
nuclei  participating  in  the  reaction  as  we  will  see  below 
in  more  detail. 

The  main  idea  of  the  QNF  procedure  is  to  approximate 
the  Hamilton  operator  iJ  by  a  simpler  Hamilton  opera¬ 
tor  obtained  from  a  power  series  expansion  of  El  which  is 
simplified  order  by  order  using  unitary  transformations. 
As  we  will  describe  in  more  detail  in  Sec.  IV  the  scaled 
Planck’s  constant,  figff,  will  play  the  role  of  a  ‘small  pa¬ 
rameter’  which  controls  the  quality  of  the  QNF  approxi¬ 
mation.  For  our  application  of  bimolecular  reactions,  the 
resulting  transformed  Hamilton  operator  truncated  at  a 
suitable  order  will  be  simpler  in  the  sense  that  it  will 


provide  an  easy,  explicit  way  to  compute  the  cumulative 
reaction  probability. 

To  define  and  implement  the  unitary  transformations 
it  is  extremely  beneficial  not  to  work  with  operators  but 
with  their  Weyl  symbols  instead.  The  Weyl  symbol  of  an 
operator  H  is  defined  as 

i7(°)(q,p;fieff)  =  y’dx(q-x/2|iJ|q-fx/2)e*P’'/''»«.  (2) 

The  superscript  (0)  is  introduced  for  reasons  that  will  be¬ 
come  clear  in  a  moment.  The  map  H  i-^  p;  hes) 

leading  to  (2)  is  also  called  the  Wigner  map.  It  is  the  in¬ 
verse  of  the  transformation  which  yields  a  Hamilton  op¬ 
erator  H  from  the  Weyl  quantization,  Op[i7],  of  a  phase 
space  function  H  (the  Weyl  map)  which,  using  Dirac  no¬ 
tation,  is  given  by 

X  y’dx|q-x/2)e-*P^/'''«(q-fx/2|. 

(3) 

Accordingly,  H^^'>  (q,  p;  hes)  in  (2)  agrees  with  the  classi¬ 
cal  Hamilton  function  H  (q,  p)  in  our  case.  The  argument 
hes  is  introduced  for  convenience  since  the  Weyl  symbol 
of  the  unitarily  transformed  Hamilton  operator  will  in 
general  explicitly  depend  on  HeS- 

We  will  now  assume  that  (q,  p;  fi.eff)  (or  equiva¬ 
lently  i7(q,  p))  has  a  (single)  equilibrium  point,  zq  = 
(qO)Po))  of  saddle-center-. .  .-center  stability  type.  By 
this  we  mean  that  the  matrix  associated  with  the  lin¬ 
earization  of  Hamilton’s  equations  about  this  equilibrium 
point  has  two  real  eigenvalues,  ±A,  of  equal  magnitude 
and  opposite  sign,  and  d  —  1  purely  imaginary  complex 
conjugate  pairs  of  eigenvalues  ±iwfc,  k  =  2, . . .  ,d.  If  the 
classical  Hamiltonian  is  of  the  form  kinetic  energy  plus 
potential  energy  then  these  type  of  equilibrium  points 
of  Hamilton’s  equations  correspond  to  index  one  saddle 
points  of  the  potential  energy.  Using  the  symbol  calculus 
the  QNF  theory  provides  a  systematic  procedure  to  ob¬ 
tain  a  local  approximation,  77qnf,  of  the  Hamiltonian  H 
in  a  phase-space  neighborhood  of  the  equilibrium  point 
Zq  in  order  to  facilitate  further  computation  of  various 
quantities,  such  as  the  CRP,  of  the  reaction  system  un¬ 
der  consideration.  In  the  following  we  summarize  the 
essential  steps  of  the  QNF  procedure. 

The  QNF  procedure  consists  of  a  sequence  of,  in  gen¬ 
eral  heS  dependent,  generalized  phase-space  coordinate 
transformations  changing  the  symbol  as 

ij(0)  ^  JjW  ^  ij(2)  ^  ff{3)  ^  ^  (4) 

The  first  of  the  transformations  (4)  shifts  the  equilibrium 
point  Zq  to  the  origin  according  to 

(z;  hes)  =  77^°^  (z  +  Zq;  Hes)  , 


(5) 
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where  z  =  (q,  p) .  Once  the  equilibrium  point  is  shifted 
to  the  origin,  the  QNF  procedure  deals  with  the  Taylor 
expansion  of  the  symbols  in  z  and  hes'. 


(z;  fieff )  =E^  +  Y.  ^ 


s=2 


with 


H. 


(n) 


(6) 


\a\  +  \l3\+2j=s  ^  10 

X  9l  ■■■(Id  Pi  ■  ■■Pd  ’ 


where  ak,Pk,j  e  Nq,  |a|  =  Y.k'^k,  |/3|  =  Lfc/^fc,  and 


H^n) 

ai,...,ad,/3i,---,/3d  J 


A  9“*^  S'®'  9'’'  /  ^ 


(8) 


At  the  next  step  of  the  transformation  sequence  one 
finds  a  symplectic  2d  x  2d  matrix  M  such  that  the  second 
order  term  of  the  symbol 

ij(2)  (z;  heff)  =  (m-1z;  H^s)  (9) 


takes  the  particularly  simple  form: 

d 

ijf  ^ (z;  heff)  =  XqiPi  +  '^Y^'^k+Pk)^  (10) 

k=2 


Section  2.3  of  Ref.  [7]  provides  an  explicit  procedure  for 
constructing  the  transformation  matrix  M . 

In  order  to  proceed  with  the  higher  order  transforma¬ 
tions  of  the  symbol  of  the  Hamiltonian  it  is  essential  to 
introduce  the  notion  of  the  Moyal  bracket.  Given  two 
symbols  A(z;  hes)  and  i?(z;  /ieff),  corresponding  to  oper¬ 
ators  A  and  B  respectively,  the  Moyal  bracket 


{A,  B}m  = 


- — A  sin 

^efT 


fteS  f  d  d 
2 


"d  ~d  \ 
dpj  dqj  j 


B 


(11) 


gives  the  Weyl  symbol  of  the  operator  i[A,  B]/heff,  where 
[•,  •]  denotes  the  commutator.  The  arrows  in  (11)  indicate 
whether  the  partial  differentiation  acts  to  the  left  (on  A) 
or  to  the  right  (on  B).  Equation  (11)  implies  that,  in 
general,  for  lieff  ^  0, 

{A,B}M  =  {A,B}^0{hls),  (12) 

where  {•,  •}  denotes  the  Poisson  bracket.  Moreover,  if  at 
least  one  of  the  functions  A,  B  is  a.  second  order  polyno¬ 
mial  in  the  variables  q,  p  then  {A,B}m  =  {A,B}.  Fi¬ 
nally,  to  simplify  further  notations  we  define  the  Moyal- 
adjoint  operator  as 


Continuing  with  the  sequence  of  transformations  of  the 
symbol  in  (4)  we  define  the  spaces 

W"  =  span  {g“G  . .  . .  ■p’^d'"  Ks  ■ 

\a\  +  \P\  +  2j  =  nj  .  (14) 

Then,  the  symbol  with  n  >  3  is  obtained  from 

means  of  the  transformation  generated  by  a 
function  lF„(z;lieff)  G  kV”) 

i7(")  =  y  ^  [Madwj'"  •  (15) 

^ '  k\ 

k=0 

The  structure  of  the  transformation  defined  by  Eq.  (15) 
implies  [7]  that  the  operators  77^"^  and  corre¬ 

sponding  respectively  (through  the  Weyl  quantization) 
to  the  symbols  and  are  related  to  one  an¬ 

other  by  means  of  the  unitary  transformation  77^”^  = 
giWnjhett  jjin-i) ^-iWn/hett  ^  whoro  Wn  is  the  operator  cor¬ 
responding  to  the  symbol  Wn-  In  terms  of  the  Taylor 
expansion  defined  in  Eqs.  (6-8)  the  transformation  intro¬ 
duced  by  Eq.  (15)  reads 


77i")=  y 


^  [Madwj''77^””^^ 
k\ 


s—k{n  —  2) 


(16) 


where  [-J  gives  the  integer  part  of  a  number,  i.e.,  the 
‘floor’-function.  Using  Eq.  (16)  one  can  show  that  the 
transformation  defined  by  Eq.  (15)  satisfies  the  following 
properties  for  n  >  3: 


77i”)  =  77i"-k  ,  for  s<n, 

(17) 

so  that,  in  particular,  =  H^\  and 

77^”)  =  77^”-^)  -  VWn  , 

(18) 

where 

7?  =  Mad^(2)  ={77f\-}. 

“2 

(19) 

Equation  (18)  is  referred  as  to  the  quantum  homological 
equation. 

We  now  specify  the  generating  function  Wn  by  requir¬ 
ing  7777^"^  =  0,  or  equivalently  77^"^  to  be  in  the  kernel 
of  the  restriction  of  T)  to  W”;  in  view  of  Eq.  (18)  this 
condition  yields 

77^n-i)  _  e  Ker  77|  w-  •  (20) 

Section  3.4.1  of  Ref.  [7]  provides  the  explicit  procedure 
of  finding  the  solution  of  Eq.  (20).  Provided  the  linear 
frequencies  a'2j  •  ■  •  j  in  (10)  are  rationally  independent, 
i.e.  m2Ui2  -I-  ...  -I-  mdUid  =  0  implies  m2  =  ...  =  md  = 
0  for  all  integers  m2,  •  ■  •  ,md,  it  follows  that  for  odd  n, 
77^"^  =  0,  and  for  even  n, 

77(")  e  span  :  |a|  +  j  =  n/2}  , 

(21) 


Mad^  :  B  1-^  MaPiaB  =  {A,  B}m  ■ 


(13) 
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where  I  =  qipi  and  Jk  =  {q\  +  Pfe)/2,  with  k  =  2, . . .  ,d, 
are  the  analogues  of  the  classical  integrals. 

Applying  the  transformation  (15),  with  the  generat¬ 
ing  function  defined  by  Eq.  (20),  for  n  =  3, . . .  ,N,  and 
truncating  the  resulting  Taylor  series  (6)  at  the  or¬ 
der  one  arrives  at  the  Weyl  symbol  H(^p  corresponding 
to  the  order  quantum  normal  form  (QNF)  of  the 
Hamiltonian  H: 

N 

i? Wp (z;  h,s)  =  Eo  +  J2  ^eff) .  (22) 

s=2 


above  and  the  recurrence  relations  (27)  and  (28).  So  the 
full  procedure  to  compute  is  algebraic  in  nature, 

and  can  be  implemented  on  a  computer.  Our  software 
for  computing  the  quantum  normal  form  as  well  as  the 
classical  normal  form  which  is  recovered  for  hes  =  0  is 
publicly  available  at  http :  / /lacms  .maths  .bris .  ac  .uk/ 
publications/ software/ index.html. 

We  stress  that  represents  an  order  approx¬ 

imation  of  the  operator  obtained  from  conjugating  the 
original  Hamiltonian  H  by  the  unitary  transformation 

fj  —  g  — iW2/fi-eff  .  .  .  q  —  IWn  /  fiell  (30) 


The  order  QNF  operator  is  then  given  by 


-^QNF  ~  -^QNF  ’ 

(23) 

where  Op  [•]  is  the  Weyl  map  defined  in  (3).  The 

Weyl  quantization  of  the  classical  integrals  /  and  Jk, 
k  =  2, .  - .  ,d,  are 

7  =  Op[/]  =  ^(qp  +  pq) , 

(24) 

Jk  =  Op[Jk]  =  ^{ql+pl),  k  =  2. 

---,d-  (25) 

Using  Eq.  (10)  and  the  linearity  of  the  Weyl 
we  get 

quantization 

d 

i/f^  =  A/  +  ^Wfcifc. 

(26) 

Since  the  higher  order  terms  in  (22)  are  polynomials  in 
I  and  Jk,  k  =  2, . . .  ,d  (see  (21)),  we  need  to  know  how 
to  quantize  powers  of  /  and  Jk-  As  shown  in  [7]  this  can 
be  accomplished  using  the  recurrence  relations 

Op  [/"+!]  =  /Op  [/"]  -  (^0  n^Op  [/"-!]  (27) 

and 

Op  [4"+i]  =  JkOp  [Jf]  +  Q)  n^Op  [jr']  (28) 


for  k  =  2, . . .  ,d.  Hence,  ^  polynomial  function 

of  the  operators  I  and  Jk- 


H, 


(N) 

QNF 


=  //g^^p(J,J2,J3,...,Jd) 

d 

=  Eq  -|-  a/  -f  ^  (  UJkJk 

fe=2 

IN/2\ 

+  E  (29) 

n—2  \ct\-\-j—n 


where  we  used  the  fact  that  the  first  two  steps  in  the  se¬ 
quence  (4)  can  also  be  implemented  using  suitable  gener¬ 
ators  Wi  and  W2  (see  [7]  for  more  details).  This  is  why 
it  is  legitimate  to  use  //qnf  instead  of  Id  in  analyzing 
such  properties  of  the  system  as  the  CRP. 

The  main  advantage  of  having  the  Hamiltonian  in 
the  form  of  a  polynomial  in  the  operators  /  and  Jk, 
k  =  2,  - .  -  ,d,  is  that  the  eigenstates  of  the  QNF  operator 
i/q^p  can  be  chosen  to  be  simultaneously  the  eigenstates 

of  the  operators  /  and  Jk,  whose  spectral  properties  are 
well  known: 

H^F\I,n2,---,nd)  =  E\I,n2,---,nd),  (31) 

where 

/|/,n2, . . .  ,nd)  =  I\I,n2,---,nd),  (32) 

Jfe|J,n2, . . .  ,nd)  =  heff{nk  +  l/2)\I,n2,---,nd)  (33) 


with  Uk  G  No  and  k  =  2,  - .  -  ,d,  and  the  energy  being 
given  by 

E  =  (^.  (n2  +  1/2), . . . ,  hesiud  +  1/2))  .  (34) 

Effectively,  the  QNF  procedure  yields  an  approximation 
of  the  original  Hamiltonian,  H ,  in  terms  of  the  operator 
i/q^p  whose  classical  counterpart  is  integrable  while  the 

classical  counterpart  of  H  is  in  general  not  integrable. 
The  approximation  is  only  valid  in  the  neighborhood  of 
the  saddle  equilibrium  point.  However,  it  is  crucial  to 
note  that  this  local  approximation  is  sufficient  to  com¬ 
pute  the  cumulative  reaction  probability  which  in  terms 
of  the  QNF  is  given  by  [7,  10] 


MiE)= 


1  +  exp  —27t 


_I{E,n2,  ---,nd) 

^eff 


n  -1 


(35) 

where  the  summation  runs  over  all  n2,  -  -  -  ,nd,  and  for 
given  energy  E  and  quantum  numbers  n2,---,nd,  the 
quantity  I  in  (35)  is  implicitly  defined  by  Eq.  (34). 


III.  COLLINEAR  HYDROGEN-  AND 
NITROGEN-EXCHANGE  REACTIONS 


The  coefficients  kn,a,j  are  systematically  obtained  by  the  gggtiou  demonstrate  the  efficiency  and  the 

QNF  procedure  to  compute  the  symbol  //q^p  desribed  capability  of  the  QNF  theory  by  applying  it  to  the  com- 
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putation  of  the  CRP  for  collinear  triatomic  reactions.  To 
this  end  we  focus  on  Hamiltonians  of  the  form 

H  =  H{qi,q2,pi,P2)  =  ^  {pj  +  pI)  +V{qi,q2) ,  (36) 

where  V{qi,q2)  gives  the  Born-Oppenheimer  potential 
energy  surface  (PES)  of  a  two-dimensional  atomic  sys¬ 
tem.  Here,  qi  and  q2  are  the  Delves  mass-scaled  coor¬ 
dinates  [12],  and  the  effective  Planck’s  constant  is  given 
by  hes  =  where  p  is  the  (dimensionless)  reduced 

mass  of  the  triatomic  system  (note  that  the  electronic 
mass  is  1  in  the  atomic  units  we  are  using). 

The  PES  is  assumed  to  possess  a  single  saddle  point 
governing  the  reaction  from  the  asymptotic  reactants  and 
products  states.  In  this  paper  we  analyze  the  following 
collinear  exchange  reactions: 

H-hH2  ^  Ha-fH,  (37) 

N-hNa  ^  Na-fN,  (38) 

where  various  isotopes  of  hydrogen  are  considered.  The 
Porter-Karplus  (PK)  PES  [13]  is  taken  to  model  the  hy¬ 
drogen  exchange  reaction  (37),  and  the  London-Eyring- 
Polanyi-Sato  (LEPS)  PES  [14]  is  adopted  for  the  nitrogen 
exchange  reaction  (38). 

We  applied  the  algorithm  presented  in  Sec.  H  to  con¬ 
struct  the  QNF  Hamiltonian  of  various  orders  for  the 
triatomic  systems  in  (37)  and  (38).  Then,  the  QNF 
Hamiltonian  was  used  to  compute  the  CRP  for  a  range 
of  reaction  energies  E  in  accordance  with  Eq.  (35).  The 
obtained  CRP-vs-energy  curves,  M{E),  were  later  com¬ 
pared  to  the  results  of  the  full  reactive  quantum  scattering 
calculations  [15,  16].  The  latter  were  performed  by  in¬ 
tegrating  the  coupled  multichannel  Schrodinger  equation 
in  hyperspherical  coordinates  [15,  16]  from  the  strong 
interaction  region  to  the  asymptotic  reactant  and  prod¬ 
uct  configurations.  The  log-derivative  matrix  method  of 
Manolopoulos  and  Gray  [17]  together  with  the  six-step 
symplectic  integrator  of  McLachlan  and  Atela  [18]  was 
used  to  integrate  the  radial  Schrodinger  equation. 

Figure  1  shows  the  CRP,  N{E)^  as  a  function  of  the 
total  energy  E  for  a  collinear  hydrogen,  ^H,  exchange 
reaction,  Eq.  (37),  on  the  PK  PES.  The  circular  points 
represent  N{E)  obtained  in  the  reactive  quantum  scat¬ 
tering  calculation,  and  can,  therefore,  be  regarded  as  the 
‘exact’  CRP  values.  The  vertical  dashed  line  shows  the 
saddle  point  energy,  Eq,  of  the  PK  PES.  The  five  solid 
colored  lines  represent  the  J^{E)  curves  corresponding  to 
different  orders,  =  2, 4, . . . ,  10,  of  the  QNF  computa¬ 
tion.  As  we  argue  in  Sec.  IV,  one  of  the  sources  of  the 
apparent  failure  of  the  QNF  method  to  reproduce  the 
correct  values  of  the  CPR  in  the  collinear  ^H  triatomic 
system  is  the  very  slow  convergence  (or  perhaps  even  di¬ 
vergence)  of  the  QNF  expansion  for  the  value  of  the  effec¬ 
tive  Planck’s  constant,  ?ieff  ~  3.07  x  10“^,  characterizing 
this  particular  reacting  system.  Another  reason  for  the 
QNF  theory  to  be  unable  to  predict  correct  CPR  values 
for  the  hydrogen  exchange  reaction  is  the  importance  of 


FIG.  1:  Cumulative  reaction  probability  as  a  function  of  the 
total  energy,  for  the  collinear  reaction  (37)  involving 

three  atoms.  The  effective  Planck’s  constant  is  fteff  ~ 
3.07  X  10“^.  The  vertical  dashed  line  shows  the  saddle  point 
energy,  Eq,  of  the  PK  PES. 

the  corner  cutting  tunneling  trajectories  [19]  in  reaction 
dynamics  of  light-atom  systems.  These  tunneling  trajec¬ 
tories  avoid  passing  through  the  immediate  neighborhood 
of  the  saddle-center-. . .  -center  equilibrium  point  in  phase 
space  and,  therefore,  their  contribution  to  the  CRP  can 
not  be  captured  by  the  QNF  theory. 


FIG.  2:  Cumulative  reaction  probability  as  a  function  of  the 
total  energy,  M{E),  for  the  collinear  reaction  (37)  with 
(tritium)  isotopes  of  hydrogen.  The  effective  Planck’s  con¬ 
stant  is  fieff  ~  1.77  X  10“^.  The  vertical  dashed  line  shows 
the  saddle  point  energy,  Eo,  of  the  PK  PES. 

Figure  2  presents  the  CRP-vs-energy  curves  obtained 
in  the  reactive  quantum  scattering  approach  (circular 
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points)  and  by  the  QNF  calculation  (colored  solid  lines) 
of  different  orders,  N  =  2, 4, . . . ,  8,  for  the  triatomic 
collinear  system  of  (tritium)  isotopes  of  hydrogen. 
The  vertical  dashed  line  shows  the  saddle  point  energy, 
Eq,  of  the  PK  PES.  The  effective  Planck’s  constant  char¬ 
acterizing  the  system  is  now  heS  «  1.77  x  10-2.  The 
convergence  of  the  QNF  ?ieff-expansion,  for  the  energies 
up  to  ~0.54  eV,  is  now  evident  from  the  figure.  How¬ 
ever,  the  QNF-predicted  CRP  values  approximate  the 
reactive  quantum  scattering  Af{E)  data  only  at  small  en¬ 
ergies.  As  in  the  case  of  the  exchange  reaction,  see 
fig.  1,  we  attribute  the  disagreement  of  the  QNF  and  reac¬ 
tive  quantum  scattering  CRP  values  to  the  non-negligible 
contributions  of  tunneling  trajectories  which  avoid  pass¬ 
ing  through  the  neighborhood  of  the  saddle. 


FIG.  3:  Cumulative  reaction  probability  as  a  function  of  the 
total  energy,  for  the  collinear  reaction  (37)  with  hy¬ 

pothetical  2*^11  isotopes  of  hydrogen.  The  effective  Planck’s 
constant  is  fieft  ~  6.9  x  10“®.  The  Af{E)  curves  obtained  with 
the  and  O***  order  QNF  are  basically  indistinguishable  for 
most  of  the  energy  range.  The  vertical  dashed  line  shows  the 
saddle  point  energy,  Eq,  of  the  PK  PES. 

Figure  3  presents  the  results  of  the  CRP  calculations 
for  a  collinear  system  of  three  hypothetical  isotopes 
of  hydrogen.  As  before,  the  circular  data  points  corre¬ 
spond  to  the  reactive  quantum  scattering  data  and  are 
treated  as  exact  CRP  values.  The  three  colored  solid 
lines  show  the  QNF  M{E)  curves  of  orders  N  =  2,4,6; 
the  M{E)  curves  obtained  with  the  4*^  and  order 
QNF  are  essentially  indistinguishable  for  most  of  the  en¬ 
ergy  range.  The  vertical  dashed  line  shows  the  saddle 
point  energy,  Eq,  of  the  PK  PES.  The  model  system  is 
characterized  by  hes  ~  6.9  x  10“^.  The  convergence  of 
the  QNF  ?ieff-expansion,  as  well  as  the  quantitative  agree¬ 
ment  of  the  QNF  predictions  and  exact  CPR  values  for 
energies  E  <  0.45  eV,  is  evident  from  the  figure. 

Comparison  of  figs.  1-3  allows  us  to  conclude  that, 
while  basically  failing  for  systems  of  light  atoms,  the 


QNF  method  of  computing  the  CPR  proves  very  effective 
for  treating  heavy-atom  reactive  systems.  On  the  con¬ 
trary,  the  full  reactive  quantum  scattering  computations 
are  only  feasible  for  reactive  systems  consisting  of  light 
atoms,  and  the  computations  rapidly  become  formidable 
as  the  atomic  mass  is  increased  [20]. 


FIG.  4:  Cumulative  reaction  probability  as  a  function  of  the 
total  energy,  J\f{E),  for  the  collinear  nitrogen  exchange  reac¬ 
tion  (38).  The  effective  Planck’s  constant  is  hes  ~  8.2  x  10“®. 
The  Af{E)  curves  obtained  with  the  4**'  and  6***  order  QNF 
are  essentially  indistinguishable  for  most  of  the  energy  range. 
The  vertical  dashed  line  shows  the  saddle  point  energy,  Eq, 
of  the  LEPS  PES. 

Finally,  in  order  to  further  illustrate  the  efficiency  of 
the  QNF  technique  for  treating  heavy-atom  systems  we 
compute  the  CRP  for  the  collinear  nitrogen  exchange  re¬ 
action  (38)  on  the  LEPS  PES.  Figure  4  compares  the 
CRP  values  obtained  in  the  reactive  quantum  scattering 
calculation  (circular  data  points)  and  those  given  by  the 
QNF  analysis  (colored  solid  lines)  of  orders  N  =  2,4,6. 
The  system  is  characterized  by  hes  ~  8.2  x  10“^.  The 
vertical  dashed  line  shows  the  saddle  point  energy,  Eq,  of 
the  LEPS  PES.  The  N{E)  curves  obtained  with  the  4*** 
and  G**'  order  QNF  are  essentially  indistinguishable  for 
most  of  the  energy  range;  this  fact  signals  the  rapid  con¬ 
vergence  of  the  QNF  ?ieff-expansion  for  the  given  value  of 
the  effective  Planck’s  constant.  The  quantitative  agree¬ 
ment  of  the  exact  and  QNF  values  of  M{E)  extends  up 
to  energies  of  ^1.5  eV. 

The  QNF  calculation  of  the  CRP  requires  signifi¬ 
cantly  less  computational  time  than  the  corresponding 
full  quantum  reactive  scattering  calculation.  For  exam¬ 
ple,  the  6*^  order  QNF  computation  of  the  nitrogen- 
exchange  CRP  curve  in  Fig.  4  took  about  10  minutes  on  a 
2.6  GHz  processor,  2  GB  RAM  computer,  while  the  cor¬ 
responding  full  quantum  reactive  scattering  computation 
took  more  than  12  hours  on  the  same  machine.  The  QNF 
approach  becomes  even  more  advantageous  for  treating 
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chemical  systems  of  atoms  heavier  than  nitrogen:  the 
expense  of  the  full  quantum  computations  rapidly  grows 
with  the  number  of  asymptotic  channels  (and,  therefore, 
with  mass)  [20],  while  the  QNF  expansion  only  becomes 
more  rapidly  convergent  making  the  corresponding  anal¬ 
ysis  computationally  cheaper. 


IV.  CONVERGENCE  OF  QNF 


While  it  is  well  known  that  for  d  =  2  degrees  of  free¬ 
dom,  the  classical  normal  form  (CNF)  converges  in  the 
neighborhood  of  saddle-center  equilibrium  points  (see, 
e.g.,  [21,  22]  )  this  is  not  clear  for  the  QNF  (for  the  hrst 
results  in  this  direction  see  [23]).  Still,  in  the  following 
we  provide  a  qualitative  discussion  of  the  convergence 
of  the  QNF  based  on  our  calculations  performed  for  the 
triatomic  collinear  reactions  of  Sec.  IIL 

The  QNF  approximates  the  Hamiltonian  of  the  reac¬ 
tion  system  in  a  phase-space  vicinity  of  the  saddle-center 
equilibrium  point.  Thus,  for  instance,  in  computing  the 
CRP  one  only  expects  this  approximation  to  render  reli¬ 
able  results  in  a  certain  energy  range  around  the  saddle 
point  energy  Eq  of  the  PES  under  consideration.  The 
energy  difference  {E  —  Eq)  may  therefore  be  considered 
as  one  small  parameter  in  the  QNF  expansion.  The  role 
of  the  other  small  parameter  is  played  by  the  effective 
Planck’s  constant,  hes-  It  is  the  convergence  of  the  QNF 
with  respect  to  this  second  small  parameter  that  we  focus 
on  in  this  section. 

We  proceed  by  considering  the  right  hand  side  of 
Eq.  (34),  i.e.,  the  QNF,  at  /  =  0,  corresponding  to  no 
‘energy’  in  the  reaction  coordinate,  and  n2  =  0,  giving 
the  zero-point  ‘vibrational  energy’  of  the  transverse  de¬ 
gree  of  freedom.  Then,  Eq.  (34)  becomes 

LAr/2j 

E  =  Eo+  Y.  CnKs-  (39) 

n—1 


For  the  case  of  the  PK  PES  the  first  five  expansion  coef¬ 
ficients  are  ci  =  0.161982,  C2  =  1.193254,  C3  =  14.90023, 
C4  =  378.7950,  and  C5  =  1227.035.  As  N  ^  00  the  radius 
of  convergence  of  the  sum  in  Eq.  (39)  is  given  by 


fc(0) 

'■eff 


■1 . 

iim  - 

n^oo  Cn+1 


(40) 


Here,  we  make  a  crude  estimate  of  by  only  consid¬ 
ering  the  first  five  expansion  coefficients  in  Eq.  (40),  i.e., 
Cn  with  n  =  1, . . . ,  5;  then,  the  radius  of  convergence  is 
given  by  ~  0.04. 

The  estimated  value  of  fifg  sheds  light  on  the  seeming 
inefficiency  of  the  QNF  theory  for  CRP  computations  in 
light  atom  reactions.  Indeed,  the  ^H  exchange  reaction, 
see  Fig.  1,  is  characterized  by  hes  =  3.07  x  10“^.  This 
value  being  close  to  fifg  signals  that  the  corresponding 


QNF  expansion  converges  very  slowly,  if  at  all,  and,  pos¬ 
sibly,  terms  of  orders  far  beyond  V  =  10  are  needed  for 
a  reliable  CRP  prediction  in  Fig.  1. 

In  the  case  of  the  ^H  exchange  reaction  the  effective 
Planck’s  constant  is  hes  =  1.77  x  10“^  and  is  thus  smaller 
than  ■  This  fact  is  in  agreement  with  the  apparent 
speed-up  of  the  convergence  of  the  CRP  values,  see  Fig.  2, 
in  comparison  with  the  ^H  case.  Finally,  the  conver¬ 
gence  is  very  fast  and  pronounced  for  the  case  of  the 
heavy  (hypothetical)  ^°H  atoms,  see  Fig.  3,  for  which 
heS  =  6.9  X  10“^  which  is  much  smaller  that  the  esti¬ 
mated  convergence  radius. 


V.  CONCLUSIONS 

In  this  paper  we  used  the  quantum  normal  form  (QNF) 
approach  to  quantum  transition  state  theory  [6,  7]  for 
computing  the  cumulative  reaction  probability  for  tri¬ 
atomic  collinear  reactions.  The  QNF  leads  to  a  realiza¬ 
tion  of  quantum  transition  state  theory  which  is  very 
much  in  the  spirit  of  (classical)  transition  state  theory. 
Similar  to  the  classical  case  where  a  recrossing  free  divid¬ 
ing  surface  can  be  constructed  from  a  classical  normal 
form  such  that  reaction  probabilities  can  be  computed 
from  the  flux  through  the  dividing  surface,  the  QNF 
can  be  viewed  to  give  quantum  reaction  probabilities  as 
the  quantum  mechanical  flux  through  the  same  (classi¬ 
cally  recrossing  free)  dividing  surface.  So  unlike  reactive 
scattering  techniques  which  involve  full,  global  quantum 
computations,  the  QNF  realization  of  quantum  transition 
state  theory  requires  only  local  information  in  the  neigh¬ 
borhood  of  the  saddle  equilibrium  point  which  governs 
the  reaction.  In  this  paper  we  demonstrated,  that  for 
heavy  atom  systems  (comprised  of  ten  or  more  nucleons) 
the  QNF  this  way  indeed  gives  a  very  efficient  method 
for  computing  cumulative  reaction  probabilities.  Here  we 
measure  ‘efficiency’  by  the  effort  for  both  implementing 
and  computing  the  QNF.  The  latter  are  both  compara¬ 
ble  to  implementing  and  computing  the  classical  normal 
form  which  lead  to  the  realization  of  classical  transition 
state  theory  (in  particular  for  multidimensional  systems). 
The  major  difference  between  the  classical  and  quantum 
case  is  that  the  QNF  computation  involves  the  Moyal 
bracket  which  slightly  more  complicated  (and  thus  com¬ 
putationally  more  expensive)  than  the  Poisson  bracket  in 
the  classical  case.  Nevertheless  the  efforts  for  implement¬ 
ing  and  computing  the  QNF  are  far  lower  than  for  the  full 
reactive  scattering  computations  to  which  we  compared 
our  results. 

We  saw,  however,  that  for  reactions  involving  light 
reactions  (such  as  the  hydrogen  exchange  reaction)  the 
QNF  gave  only  very  poor  results.  We  attributed  the  fail¬ 
ure  of  the  QNF  computation  in  these  cases  to  the  pres¬ 
ence  of  corner  cutting  tunneling  trajectories  which  are 
not  captured  by  the  QNF.  This  way  the  QNF  and  reac¬ 
tive  scattering  methods  can  be  viewed  as  complementary 
methods  where  the  latter  gives  very  good  results  for  light 


atom  systems  and  the  former  displays  its  full  power  espe¬ 
cially  for  heavy  atom  systems  for  which  reactive  scatter¬ 
ing  approaches  become  very  difficult  or  even  unfeasible 
due  to  the  growing  number  of  reactive  channels  that  have 
to  be  taken  into  account  [20]. 

We  note  that  also  other  approximation  techniques  such 
as  the  initial  value  representation  (IVR)  [24]  have  been 
shown  to  be  fruitful  for  reaction  probability  analysis  of 
collinear  triatomic  reactions  [25].  However,  in  order  to 
properly  account  for  interference  effects  the  IVR  method 
requires  propagation  of  a  huge  number  of  classical  trajec¬ 
tories  and,  therefore,  can  pose  difficulties  for  application 
to  high-dimensional  atomic  systems  whereas  the  difficul¬ 
ties  in  computing  the  QNF  do  not  grow  so  rapidly  with 
the  number  of  degrees  of  freedom.  In  fact  it  would  be 
very  interesting  to  make  a  detailed  comparison  between 
the  QNF  and  the  IVR  approach. 

Another  benefit  of  the  QNF  approach  to  compute 
cumulative  reaction  probabilities  lies  in  the  fact  that 
it  involves  only  little  (local)  information  of  the  Born- 
Oppenheimer  PES;  namely  the  Taylor  expansion  of  the 
PES  about  the  saddle  equilibrium  point  governing  the 
reaction.  In  fact  we  saw  that  highly  accurate  results  over 


quite  a  broad  energy  range  can  already  be  obtained  from 
the  4*^  or  6*^  Taylor  expansion  which  enters  the  QNF 
of  the  same  order.  This  is  especially  useful  for  systems 
for  which  the  computation  of  the  global  PES  required  in 
other  methods  is  very  difficult. 
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